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Abstract 

Electromigration phenomena in metallic lines are studied by using a biased 
resistor network model. The void formation induced by the electron wind 
is simulated by a stochastic process of resistor breaking, while the growth 
of mechanical stress inside the line is described by an antagonist process of 
recovery of the broken resistors. The model accounts for the existence of 
temperature gradients due to current crowding and Joule heating. Alloying 
effects are also accounted for. Monte Carlo simulations allow the study within 
a unified theoretical framework of a variety of relevant features related to the 
electromigration. The predictions of the model are in excellent agreement 
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with the experiments and in particular with the degradation towards electrical 
breakdown of stressed Al-Cu thin metallic lines. Detailed investigations refer 
to the damage pattern, the distribution of the times to failure (TTFs), the 
generalized Black's law, the time evolution of the resistance, including the 
early-stage change due to alloying effects and the electromigration saturation 
appearing at low current densities or for short line lengths. The dependence of 
the TTFs on the length and width of the metallic line is also well reproduced. 
Finally, the model successfully describes the resistance noise properties under 
steady state conditions. 



PACS numbers: 66.30.Qa; 85.40.Qx; 85.40.Ls; 64.60.Ak 
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I. INTRODUCTION 



The phenomenon of electromigration (EM) is typical of metaUic conductors and consists 
in a non-steady atomic transport driven by electronic currents of high density^'^. The non- 
steady atomic transport gives rise to the formation and growth of voids and hillocks in 
different regions of the conductor (EM damage) This damage cumulates progressively 
and, when the current is applied for a sufficiently long time, a void grows enough to break 
completely the metallic line implying an irreversible failure process. The time required for 
this process defines the time to failure (TTF) of the metallic hne^'^ (though alternative 
failure criteria can be found in the literature^^^). The importance of EM is largely due 
to the fact that it is the most common mechanism of failure of the metallic interconnects 
present in any electronic device^"^. As a consequence, a huge number of experimental^"^^ 
and theoretical^'''"^^ studies have been and are yet devoted to the subject, especially in the 
context of modern nanoelectronics. 

The central issue in EM degradation phenomena is the determination of the TTF and 
its statistical properties'"^. TTFs are measured under very high stress conditions (currents 
and temperature much higher than those corresponding to the usual operating conditions 
of the devices) in the so called accelerated tests''^. The extrapolation to normal operating 
conditions is generally performed on the basis of three frequently adopted assumptions^'^. 
First, TTFs are taken to follow a lognormal distribution (which means that the logarithms of 
TTFs are normally distributed). This distribution is then characterized by two parameters: 
the median time to failure, t^o, and the shape factor, s, where t^o is the time corresponding 
to the failure of 50% of the lines in the statistical sample and s is the lognormal root-mean- 
square deviation'. This distribution has been observed in many EM failure tests'"^'^'^''^ 
and recently, new testing techniques'^ allowing the analysis of very large statistical samples, 
have shown that the TTF distribution follows a perfect lognormal behavior down to four 
shape factors'^. In spite of this evidence, no satisfactory eplanation has been given until 
now for the lognormality of the TTF distribution''^''^. 
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The second assumption concerns the independence of the distribution shape factor of the 
stress conditions. Actually, s was found to be independent of temperature in all the range 
of values usually considered in accelerated tests^^^'^'®'^^. On the contrary, a broadening of 
the distribution has been observed at the lowest current densities used in these tests^'^'^. 
This broadening of the distribution at low stress conditions has crucial implications on the 
evaluation of the minimum time to failure, i.e. the time corresponding to the first failure of 
a line of the family^'^. We will discuss this point in Sec. III.B. together with the results of 
our simulations. 

The third assumption concerns the validity of the following empirical law, known as 
Black's law^'^'^^, relating to the current density, j, and temperature, T: 

i5o = Cj-"exp[^] (1) 

where C is a fitting amphtude, n the so called current exponent, E the EM activation 
energy and ks the Boltzmann constant. The validity of this law is confirmed by many 
experiments^^^. However, the simple power law behavior described by Eq. (1), generally 
holds in a limited range of intermediate current densities^^^. At high current densities, 
another power law with a different (higher) exponent is frequently observed and usually 
attributed to Joule heating effects^^^'^^, while at low current densities t^o deviates from 
a power law, showing a tendency to diverge^'^. Furthermore, even in the intermediate 
range of j values, there is no agreement in the literature about the value of the current 
exponent^^^'^'''^^. We will come back on this subject in Sec. lll.B, discussing our results. 

The geometry of the line plays a crucial role on the EM damage^"^'^'^^. In fact, the deple- 
tion and accumulation of mass in different regions of the film, under the driving force exerted 
by the electronic current, determines the growth of mechanical stress gradients. The last 
ones give rise to an atomic back- flow which contrasts the EM process^'^'^^. As the strength 
of these gradients depends strongly on the line geometry^'^'^'^^, geometrical parameters have 
a fundamental role in the occurrence of failures. In particular, the competition between the 
driving force of the electronic current and the action of mechanical stress gradients, results 



in the existence of a current density threshold below which EM is stopped^^. This condition 
was expressed by Blech and Herring^^ in the following equation which relates the threshold 
current density with the length of the line: 

where Qa is the atomic volume, Act is the maximum value of the mechanical stress difference 
between the line terminals that a line of length L can bear, p is the resistivity of the line, 
e the electronic charge and Z* an average effective charge^'^'^^. Equation (2), known as 
Blech's law^'', defines the so called threshold product, {jL)c. We will discuss the role of the 
line geometry on the EM failure process in Sec. III.C. 

Another fundamental ingredient in the understanding of the EM damage of interconnects 
is represented by the granular structure of the materials employed, Al, Cu, Ag, Al alloys, 
etc^~^. Furthermore, it must be noted that a high degree of disorder is usually present in alloy 
films (typically Al-Cu, Al-Si) due to alloying effects^'^^'-*^^'^^"^" and to thermal gradients^'^. 

A large number of models has been proposed for the study of EM^^~^^. Many of them are 
microscopic models which address the problem of identifying the mechanisms responsible for 
the degradation process in terms of the pecuharities of the material considered^^'^^'^^'^^'^^. 
Then, by using appropriate kinetic equations, some specific features of the damage process 
are determined and compared with experiments^^'^^'^^'^^'^^. If the peculiarities of the material 
are sufficiently well accounted for, the predictivity of these models is very high^^'^^'^^ and 
therefore they can be very useful for applicative purposes. However, the approach used by 
many of these models intrinsically limits their predictivity to some specific features of the 
EM damage. 

Another class of EM models has also been developed in the literature^°'^^, based on 
a "coarse grain" random resistor network (RN) approach^^^^^. Actually, the use of this 
kind of models is particularly appropriate in consequence of the granular structure of the 
materials used for the interconnects. Indeed, it has been observed that the atomic transport 
through grain boundaries and interfaces (transport channels) far exceeds that through the 



bulk of the grains^'^. Therefore, it is generally possible to neglect mass transport everywhere 
except within these channels and to describe the film as an interconnected network of atomic 
conducting paths^. 

Bradley et al.^° were the first to propose and apply to the study of EM a kinetic version 
of the random fuse modeP^. The dynamic fuse model introduced by these authors^'^ adopts 
a failure criterion for the elementary resistor of the network suitable for the description of 
the EM process. The predictions of this model concerning the damage pattern and TTFs 
(maximum and minimum TTF, relationship of the median time to failure with the current, 
the temperature and the line geometry) shows the effectiveness and the potentiality of the 
resistor network approach^°. However, the model of Bradley et al. cannot describe many 
other features of EM. In fact this model, though giving a good description of the driving 
force of the electronic current, does not take into account the antagonist action exerted by 
mechanical stress gradients. On the other hand, the competition between these two effects is 
essential for giving rise to a threshold current for EM^'^. Therefore, all the phenomenology 
related to the Blech's law^'^'^'^^ and saturation effects^ cannot be accounted for within the 
Bradley's model. Moreover, this model completely neglects Joule heating effects which are 
present in high stress condition^'^. 

Here, we illustrate a theoretical approach to EM which aims at studying the different 
features associated with this phenomenon within a unified theoretical framework. Similarly 
to the approach of Bradley et al. , our study is performed by renouncing to provide a descrip- 
tion of the kinetics at an atomistic level and by adopting the RN approach, thus focusing on 
the correlations established by the electronic current among the different components of the 
system (grains, clusters of grains, interfaces, etc. i.e. atomic transport channels). However, 
in contrast to the dynamic fuse modeP'^, we use the biased percolation model^^'^^, which 
adopts a probabilistic failure criterion for the elementary resistor of the network. The EM 
damage is described in terms of competition between two biased stochastic processes taking 
place in a resistor network^^. Then, by means of Monte Carlo (MC) simulations, we are 
able to study a variety of relevant features of EM degradation. Early stage results has been 
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presented in Refs. [ 28, 38]. In this article, we present a comprehensive study which includes 
many fundamental new features. 

As a first we will show results concerning the damage pattern, the resistance evolution and 
alloying effects. Then, large attention will be devoted to the behavior of TTFs. We will show 
that the model correctly predicts a lognormal distribution for them, perfectly superimposing 
with the experimental one. The dependence of the parameters of the TTF distribution on 
temperature and current has also been investigated. The Black's law^'^'^^ has been recovered 
not only for low but also for high current values and the current exponents are in good 
agreement with the experimental results in both cases. By using a rectangular geometry, 
the dependence of TTFs on the length and width of the metaUic hne has been investigated, 
so that, for a fixed width, the Blech's law^'^'^^ has been tested and a general expression has 
been obtained for the dependence of TTFs on both, length and width. Finally, but not 
of minor interest, we have considered resistance saturation effects^ and the properties of 
resistance fluctuations, by focusing on the non-Gaussianity of the distribution and on their 
power spectrum. Thus, our approach is able to account for a phenomenological scenario 
much wider than that considered by all existing EM models. 

Finally, we emphasize the fact that the interest in EM phenomena is not limited to strictly 
applicative and practical purposes. Indeed, the understanding of non-equilibrium and failure 
phenomena in disordered systems represents a fundamental topic which has attracted a large 
attention in the recent hterature^^"^^'^^"^^. In this respect, it must be noted that EM, which 
occurs in granular materials, in presence of a significant disorder, driven by an external bias 
and contrasted by growth of mechanical stress gradients, exhibits practically all the main 
ingredients to represents a paradigmatic example of failure process in a disordered system. 

The paper is organized as follows. Sec. II briefly surveys the theoretical model and 
define the parameters of interest. Sec. Ill presents the results in connection with: (A) 
resistance evolution, (B) stress conditions, (C) geometrical effects, (D) resistance saturation 
and fluctuations. Sec. IV draws the main conclusions. 
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II. MODEL 



We describe a thin metallic line of length L, width W and thickness th as a two- 
dimensional resistor network of rectangular shape and square-lattice structure. The net- 
work of resistance R is made by Nl and Nw resistors in the length and width directions 
respectively. The external bias, represented by a constant voltage y or a constant current 
/, is apphed to the RN through electrical contacts realized by perfectly conducting bars at 
the left and right hand sides of the network. Thus, the total number of network resistors 
(excluding the contacts) is: Nfot — 2NlNw + Nl — Nw Each resistor can be associated 
with a single grain, a small cluster of grains or an interfacial path. By denoting with d the 
average size of the grains, of the grain clusters or of the interfacial path, the values of Nl 
and Nw can be related to the ratios L/d and W/d, respectively. The network lies on an 
insulating substrate at temperature Tq, acting as a thermal bath, and it is made by three 
kinds of resistors: i) regular resistors, ii) impurity resistors, iii) broken resistors. The regular 
resistors are associated with grains of "normal" resistivity (void free). The resistance of 
these resistors depends linearly on temperature, according to the expression: 

rregATn) = ^re/[l + "(^n " ^re/)] (3) 

Where n is the resistor label, a the temperature coefficient of the resistance (TCR), T„ 
the local temperature, T^e/ and r^e/ the reference values for the TCR. When the Joule 
heating is neghgible T„ = To and the regular resistors are all equal to ro = rreg{To). The 
impurity resistors of resistance r^mp < ''"o are associated with the formation/dissolution of 
CuAl2 precipitates (low-resistivity cluster) during the stress conditions. Thus, they account 
for the variation in the alloy composition during the EM test due to the electron wind 
and/or to thermal effects (alloying effects) Finally, the broken resistors correspond 

to the presence of microvoids at the grain boundary and, possibly, inside the grains. These 
broken resistors of resistance rop (OP stays for open circuit) are thus associated with very 
high resistivity regions inside the line. Here, we have taken rop — lO^ro. The existence of 
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temperature gradients due to current crowding and Joule heating effects is accounted for by 
taking the local temperature of the n-th resistor (regular, impurity or broken) of resistance 
r„ given by the following expression^^: 



Tr, = To + A [rnil + — Yl {^m,n^m,n " ^nil)] , (4) 

where A is the thermal resistance of the single resistor, Nneig is the number of first neighbors 
of the n-th resistor, i„ the current flowing in it and im,n is the current flowing in the m-th 
neighbor. The value B — 3/4 is chosen to provide uniform heating of the perfect network, 
i.e. made by identical resistors. Equation (4) is the Fourier equation written by taking 
the simplifying assumption of instantaneous thermalization of the resistor, i.e. by taking a 
stationary regime and neglecting time dependent effects in the heat diffusion^^. Under this 
hypothesis, the diffusive term must balance the term related to the power supplied by the 
external current and generated on the resistor or transmitted to it by mutual interactions. A 
simplified form of Eq. (4) , not including thermal exchanges with neighbor resistors has been 
firstly proposed by Gingl et al.^^. For a perfect or nearly perfect network of resistance Rq, 
when a mean-field approach is meaningful, the average temperature increase is thus^^: 

AT = ARoiyNtot = eRol\ (5) 

where 9 — A/Nfot is the structure thermal resistance^'^. 

The EM damage, consisting in the formation of microvoids under the action of the 
electron wind, is simulated by a stochastic process of resistor breaking. In other terms, we 
consider that the transformation of the n-th resistor (regular or impurity) into a broken one, 
T^n fop-i can occur with probability Wop,n- By adopting the biased percolation modeP^'^^, 
we have taken the following expression for Wop,n- 

WoP,n = ^W[-Eop/kBTn]-, (6) 

where Eqp is a characteristic activation energy. In fact. Equation (6) coupled with Eq. (4), 
implies that the void formation process is a biased percolation^^'^^. This means a probability 



of breaking a resistor (generation of microvoids) higher for resistors crossed by high current 
values^^'^^. Thus, the breaking probabihty is nonuniform for the different resistors in the 
network and it changes with time. We note that Wop,n depends on the current distribution, 
which in turn depends on the network configuration^^'^^. As the last one results from a 
progressive accumulation of damage, the history of the network is partly accounted for by 
the biased percolation model. 

The effect of the atomic back-flow, which contrasts the EM, is simulated by introducing 
a recovery process consisting in a stochastic heahng of the broken resistors^^'^^'^°. Therefore, 
the transformation rop ^ Treg,n is allowed with probability: 

WR,r^ = ^M-ER|kBT^i (7) 

where Er is a recovery activation energy. Furthermore, the variation of the line composition 
due to alloying effects is described by allowing the stochastic transitions^^: rreg,n fimp 
and Timp rreg,n, occurriug with probabilities Wm^n = (^wi-Em/kBTn] and Wm^n = 
exp[—Eiji/kBTn] respectively, where Eri and Ej^ are two characteristic activation energies. 
By adopting this description of alloying effects, we are limiting ourself to consider the change 
in the elementary resistances due to a variation inside any grain or cluster of grains (single 
resistor) of the number of Cu atoms dispersed in the matrix or present in small CuAl2 
precipitate^^. 

The initial configuration of the RN can contain some initial concentrations of broken and 
impurity resistors, Pini and p^^^, respectively. Alternatively this initial configuration can be 
chosen as the perfect one: pini — and p^^l — . The network evolution is obtained by 
Monte Carlo simulations which are carried out according to the following iterative procedure. 

(i) Starting from the initial network, we calculate {in} and the network resistance R by 
solving Kirchhoff's loop equations. Moreover, we calculate {Tn} by using Eq. (4). 

(ii) OP and rimp are generated with the corresponding probabilities Wqp and Wm and 
the remaining Vreg are changed according to {Tn}. Then {in} and {Tn} are recalculated. 

(iii) OP and rimp are recovered with probabilities Wr and Wjr respectively, and the 
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resistances r^eg are changed again according to 
(iv) {in}, {Tn} and R are recalculated. 

This procedure is iterated from (ii), thus the loop (ii)-(iv) corresponds to an iteration 
step, which is associated with a unit time step on an arbitrary time scale to be calibrated by 
comparison with experiments. Depending on the parameter values, which are related to the 
physical properties of the line and to the external conditions, the two following possibilities 
can be achieved during the iteration procedure: irreversible failure or steady-state evolution. 
In the first case, at least one percolating cluster of broken resistor (i.e. a cluster connecting 
the top and the bottom of the network) is formed and thus the resistance R diverges^^. In the 
second case, the network resistance fluctuates around an average value < R > (saturation 
value) The average over the statistical ensemble (different realizations of failure of 
networks with the same parameters and in the same external conditions) of the values of 
the fraction of broken resistors corresponding to the appearance of at least one percolating 
cluster, is called percolation threshold and it is denoted as Pc^^- 

To check the model we have considered EM tests performed with a standard median 
time to failure technique^'^ on Al-0.5%Cu hues. The tests have been carried out at different 
currents and temperatures by adopting a 2-metal level configuration with tungsten vias^'^ 
and by using a 20% relative resistance variation as failure criterion. The lines used in the 
tests were 3000 iim long, 0.45 iim wide and 0.8 //m thick. The last thermal treatment 
undergone by these lines occurred during fabrication and it consisted in a high temperature 
annealing followed by a rapid cooling. This treatment left a non-equilibrium concentration 
of Cu dissolved into the Al matrix. Therefore, in the early stage of the EM test, the heating 
associated with the stress conditions gives rise to a formation of CuAl2 precipitates. The 
low resistivity of these clusters and, mainly, the reduction of the internal disorder, i.e. the 
reduction of the number of scattering centers of Cu in the solid solution, cause an initial 
decrease of the hue resistance^^. Here, we consider the data obtained at T = 492 K, shown in 
Ref. 28, and the data obtained at T = 467 K, reported in the following section. In both cases 
the stress current density was j — S MA/cm^. The resistance of the lines at the reference 
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temperature T^e/ = 273 K was: = 269 Q (averaged over a family of 40 samples), while 
the TCR was 3.6 x IQ-^ K'^. 

The values of the parameters used in the simulations have been chosen as follows. We 
have taken the values corresponding to the actual stress conditions and to the physical 
parameters of the metallic line whenever possible, i.e. whenever it was possible to make 
a direct correspondence with the model parameters and the line properties and when this 
choice was not too heavy computationally. The remaining parameters have been chosen to fit 
the experimental results and/or to reduce the computational effort. Concerning this point, 
we notice that the present approach allows for a direct simulation of lines with ratio L/W up 
to 150. To describe the resistance evolution of lines characterized by higher values of this 
ratio, as the lines tested in the experiments shown in the next section (where L/W — 6667), 
we have adopted the following further approximation. The network is taken to represent 
the region of dominant void growth inside a longer line, i.e. the region responsible for the 
resistance variation of the hne. We can take the length of this region given by L/F, where 
F is an integer number. Thus, in the initial conditions: Ro,iine = FRq. On the other 
hand, according to the above assumption we have: ARune — ^R- Therefore, the relative 
resistance variation of the whole hne can be expressed as: ARune/ Ro,iine — {^/ F){A.R/ Rq). 
We underline that this approximation is used only to check the model by a direct comparison 
of the measured resistance evolutions of long lines (Fig. 2) and the evolutions calculated by 
the present model (Fig. 3). All the other results concern short hues and do not make use of 
the above approximation. 

Thus, except when differently specified, we have used the following values of the param- 
eters: Nw = 12, Nl = 400, F = 1, To = 492 K, 7 = JWth = 10.8 mA (which corresponds 
to the values of j, W and tth used in the EM tests cited above), a = 3.6 x 10^'^ K~^, 
Tref — 273 K, Tref — 0.048 fl, Timp — 0.016 fl. FoT the long lines used in the EM tests, when 
F — 200, this value of r^e/ provides the correct value for R^ef, reported above. Moreover, 
we have taken A = 2.7 x 10^ K/W. According to Eq. (5), this value of A provides an initial 
heating of the network of 8.3 K, comparable with that estimated in the experiments of Ref. 
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28. The initial network configuration corresponds to pini — and p-^^ = 0. Furthermore, 
we have taken: Eqp = 0.41 eV and Er = 0.35 eV, as reasonable values for the activation 
energies Eqp and Er. We notice that the value of Eqp controls the range of the time scale, 
nevertheless this range is in any case arbitrary within our model. Thus, the value of Eqp 
can be considered as a free parameter which can be chosen with the purpose of saving com- 
putational time. More crucial is the choice of E'ij, whose value sets the importance of the 
recovery process^^, i.e. the strength of the atomic back-flow due to the mechanical stress. 
Consequently with this choice of Eqp and Er, the values of Erj and Eir have been taken 
sufficiently small to account for the separation of the temporal scales of the void formation 
and of the alloying processes, observed in the experiments. Here, we have taken: Eri — 0.22 
eV and Eir = 0.17 eV. 

III. RESULTS 

A. Resistance evolution 

A typical resistance evolution and the corresponding damage pattern near the final failure 
are reported in Figs. 1(a) and 1(b), respectively. In this case Nl = 48 while all other 
parameters take the values specified at the end of the previous section. In Fig. 1(a) we 
observe a resistance drop at the early stages of the evolution due to the generation of 
impurity resistors. This process rreg,n ~^ ^imp simulates the variation of composition of the 
line associated with the initial precipitation of part of the Cu dissolved in the Al matrix, 
as discussed in Sec. II. On the other hand, this process together with the antagonist one, 
^imp ^reg,n) takcs placc ou a timc scale that is much shorter than that associated with 
the generation of broken resistors. Therefore, after a given amount of time (relaxation time 
of the alloying process) , the concentration of the impurity resistors reaches its steady state 
value^^ corresponding to the temperature Tq + AT. On a longer time scale, the fraction of 
broken resistors increases and the network becomes more and more unstable. This implies 
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an increase of both the resistance value and the resistance fluctuations, as shown in Fig. 
1(a). Finafly, at a given time (time to failure) the fraction of broken resistors reaches the 
percolation threshold and R{t) diverges. Figure 1(b) reports the damage pattern just before 
the flnal failure. Precisely, this flgure shows the temperature distribution inside the network: 
the broken resistors are the missing ones while the different gray levels, from black (cold) to 
white (hot), correspond to different T„ values ranging from the substrate temperature up to 
700 K, with a temperature step of 10 K. The damage pattern mainly consists of a channel of 
broken resistors elongated in the direction perpendicular to the current flow, a characteristic 
feature of the biased percolation^^'^''. This simulated damage pattern reproduces well the 
experimental pattern observed by scanning electron or x-ray microscopy in metallic lines 
which are failed due to EM^'^'^^'^^. 

Figure 2 shows the resistance evolutions of seven Al-0.5%Cu lines measured in the EM 
tests performed at T = 467 K, as described in Sec. II. These lines were stressed by a 
current density j — S MA/cm^ which corresponds to / = 10.8 mA. Figure 3 reports the 
resistance evolutions obtained by simulations. Here, different curves correspond to different 
reahzations of failure. In this case, we have taken: Nw — 12, Nl — 400, F — 200, Tq = 467 
K, Tref = 0.044 Q, Timp = 0.006 Q, Pi^i = (2.5 ± 0.2) X 10~^ (the broken resistors in 
the initial network configuration are supposed uniformly distributed), while the remaining 
parameters have the same values specified at the end of Sec. II. The time scale in Fig. 3 has 
been calibrated according to the following procedure. The statistical sample tested in the 
experiments was composed by thirteen lines and the resulting median time to failure, ^50, 
was: t^o^exp ~ 1-3 x 10^ s. A sample of thirteen simulated failure realizations was considered 
and the corresponding t^o^sim was calculated in units of iteration steps. From these two 
values, we obtained the value At = 185 sec for the time interval to be associated with each 
iterative step. The comparison between Figs. 2 and 3 shows that the calculated evolution 
of the resistance well reproduces the main features of the observed evolution. We note that 
a small discrepancy between measured and simulated evolutions appears just before the 
failure, where the abrupt increase of the resistance shown by the simulated curves, contrasts 
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with a pre-failure increase of the resistance generally present in the experimental ones. This 
discrepancy, is partly due to the factor 1/F relating the relative resistance variation of 
the network with that of the metallic line. In fact, the simulated evolutions obtained for 
short hues, see Fig. 1(a), when the full metallic line can be directly simulated and the 
approximation adopted for long lines can be avoided (i.e. F = 1) display also a pre-failure 
region. 

The agreement is further confirmed by the comparison between the distributions of the 
measured and calculated TTFs reported in Fig. 4. This figure shows on a lognormal plot^^ 
the cumulative distribution function (CDF) of the failure probability as a function of the 
times to failure obtained from the EM tests (full circles) and simulations (open circles) 
considered above. The agreement between experiments and simulations is excellent and the 
shape factor of the two distributions is s = 0.16 in both cases. Thus, the direct comparison 
of the results of the model and of the EM tests shown in Figs. 2, 3 and 4, vahdate the 
present computational approach for the study of EM failures. 

As anticipated in Sec. I, two central problems encountered in the study of EM phe- 
nomenon concern the role played by the stress conditions and the line geometry on the 
damage process^'^'^^. Therefore, with the purpose of further checking the predictivity of the 
model and of extracting new information from it, we have calculated the effect on TTFs of 
temperature, stress current, length and width of the lines. To contain the computational 
effort and to avoid the approximation used for long lines, the study has been limited to short 
or moderately short lines. Thus, in the following we will discuss the results of simulations 
carried out by taking F = 1, by varying Tq, /, Nl, Nw, while keeping all the remaining 
parameters to the values specified at the end of Sec. II. 

B. Stress conditions: temperature and current effect 

We start by considering the effect of temperature on the times to failure of 20 networks of 
sizes 12 X 400 stressed by a current of / = 10.8 mA. Accordingly, we analyse the dependence 
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on temperature of ^50 and s, the two parameters which determine a lognormal distribution. 
We consider fourteen values of Tq ranging from 400 K to 800 K, a considerably wider range 
with respect to standard accelerated tests^'^. Figure 5 shows the cumulative distribution 
functions of the failure probabihty calculated for To = 800 K (triangles left), Tq — 650 K 
(open squares), Tq — 467 K (open circles) and Tq = 400 K (triangles down), while the 
solid lines fit the CDFs with lognormal distributions. These four Tq values are selected as 
representative for the behavior of s. Indeed, we have found that s is nearly independent of 
temperature in a wide range of intermediate temperature values, while it increases signifi- 
cantly at low temperatures and decreases at the highest temperatures considered here. The 
broadening of the distribution at low temperatures and its narrowing at high temperatures 
witness the different importance of the network microgeometry in the two extreme stress 
conditions of 400 and 800 K, respectively. Indeed, the network microgeometry, resulting 
from the stochasticity of the defectiveness, give rise to a kind of network individuality. At 
very high stress, the differences in the network microgeometry loose their importance. The 
contrary occurs at low stress, where this diversity becomes of importance in determining the 
actual TTF. We underline that the broadening of the TTF distribution at low temperatures 
has important implications in the interpretation and use of the results of accelerated EM 
tests. In fact, when evaluating the reliability of a family of lines under standard operating 
conditions (usually close to room temperature and relatively low current density) , it is cru- 
cial to estimate not only but also the minimum time to failure^'^. The determination of 
these two quantities is usually obtained from accelerated tests performed at high tempera- 
tures on a statistically significant, but in any case small, sample of the entire family. Then, 
the estimate of the minimum time to failure of the family is obtained by an extrapolation 
of the CDF in the region of low failure probability^'^. Such an estimate is very sensitive 
to a possible broadening of the TTFs distribution at the operation temperature. For this 
reason it is crucial to estimate and to take into account the dependence of s on temperature. 
We remark that the increase of s at low temperatures is a source of the following apparent 
paradox®: the minimum time to failure at low temperatures can be shorter than the mini- 
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mum time to failure at high temperatures. We will face a similar paradox by discussing the 
dependence of s on the current^. Solutions to this apparent paradox has been proposed in 
the literature®' based on the necessity of testing large samples and on the introduction 
of a three parameters lognormal distribution, where the third parameter is a characteristic 
incubation time®. 

The analysis of the temperature dependence of the TTFs is completed by Fig. 6 which 
displays on linear-log scale the calculated values of ^50 as a function of the inverse of the 
substrate temperature. Here, the dashed line is the best fit obtained with the function 
Z exp[E /UbTq], where Z is & fitting amplitude. Thus, the calculated values of ^50 perfectly 
follow, within the numerical uncertainty, the Black's law^'^'^^, discussed in Sec. I. The value 
of E extracted from the fit is £J = 0.41 eV, thus E — Eqp-, and we can identify the activation 
energy of the resistor breaking process with the EM activation energy. 

To investigate the effect of current on the failure process, we have calculated the TTFs 
of 12 X 400 networks stressed at Tq = 492 K by different (thirteen) current values in the 
range 5.0 -^ 60. mA. For each current value s and tso have been determined by considering 
twenty realizations of failure. Figure 7 shows the CDFs of the failure probability versus 
failure time, calculated for I — 7.5 mA (triangles down), I — 30. mA (open squares), 
I — 60. mA (open circles). The solid lines fit the CDFs with lognormal distribution. We 
have found that the shape factor of the distribution exhibits a minimum at intermediate 
values of /. Thus, we can identify two regions of current values: a moderate current (m.c.) 
region, where s decreases at increasing current, and a high current (h.c.) region, where s 
increases at increasing current. Such a broadening of the TTFs distribution at low currents 
has been actually observed in several EM tests^'®. Its imphcations concerning the evaluation 
of reliability of metallic lines are similar to those previously mentioned in connection with 
the effect of temperature. A detailed discussion of these problems can be found in Ref. 
8. Here, we underline that the non-monotonic behavior of s versus current contrasts the 
general monotonic behavior found versus temperature. This non-monotonic behavior of s 
can be understood by considering also the dependence of ^50 on the current that is reported 
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in Fig. 8. More precisely, we show in the inset of this figure a log-log plot of the calculated 
values of t^o versus /. We can see that exhibits two power-law regions. A first one 
is in the moderate current region and is characterized by a current exponent n = 2.1, 
in good agreement with the Black's law^^'^^. A second power-law with a higher exponent, 
n = 5.7, is found in the high current region. We notice that similar behaviors at high current 
densities have actually been observed in EM measurements and are frequently reported in 
the hterature^"^. Furthermore, the same inset of Fig. 8 shows that ^50 drastically increases 
at the lowest currents. Here, for / — 4.0 mA the networks are found to remain stable over 
more than 5 x 10^ iterations. This strong increase of ^50 is associated with the existence of 
a threshold current, Ib, below which a steady state condition is achieved, manifesting itself 
in a saturation of the network resistance. Accordingly, for I < Is the electrical breakdown 
no longer occurs. For this reason, in Fig. 8 the region corresponding to I < Ib is evidenced 
in gray. Previous investigations of the general properties of the model, reported in Refs. 49, 
have shown the existence of this threshold. Some of the properties of the steady state of the 
network will be discussed later in connection with the results concerning resistance noise and 
we refer the reader to Refs. 49 for a deeper analysis of these properties. We emphasize that 
this sharp increase of ^50 at low currents has also been observed in EM tests^'^. Therefore, the 
dependence of ^50 on the current obtained by simulations agrees with the behavior measured 
over the full range of current values. To complete the analysis. Fig. 8 reports in a log-log 
plot the median time to failure versus the difference I — Ib (full squares) . By taking for Ib 
the value given above, Ib = 4.0 mA, we have found that t^o scales as: 

t5o-(/-/s)-"« (8) 

This expression, firstly proposed by Filippi et al.^, can be considered as a generahzation 
of the Black's law and Ug can be called as generalized current exponent, to distinguish it 
from the current exponent n of the conventional Black's law, Eq. (1). Figure 8 displays two 
current regions, each characterized by a given value of Ug. These two regions correspond to 
the different current dependence of the distribution shape factor, reported in Fig. 7. The 
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common behavior of t^o and s with current originates from the change with the bias of both 
the damage pattern and the magnitude of Joule heating and can be understood as follows. 

First, let us consider a simpler system: a network in which two stochastic processes of 
resistor breaking and recovery occur with uniform probabilities: Wdo — exp[— £'op//cb3o] 
and Wro — exp[— _Er/A:bTo]. A similar network, subjected to random percolation^^'^^ de- 
scribes well the instability of very thin metallic films due to agglomeration phenomena^^. 
The stability of this network has been studied in Ref. 52, where the failure condition and 
the expression for the average time to failure^^ (ATTF) have been derived in the limit of 
networks of infinite size. Here, to point out the dependence of the percolation threshold on 
the system size, it is convenient to write the failure condition in the following form: 

similarly, the average time to failure can be written as: 

H(1-Wdo){1-Wro)] ^ > 



with: 



Pc 



1 + 1 (11) 



where g < 1 for failing networks, according to Eq. (9). For this simple system it is quite 
easy to determine the role of the size and of the geometry of the network on the value of the 
percolation threshold. In fact, an ideally infinite network [N^ — > oo and Ny/ oo), with 
percolation threshold Pc,oo, would have a zero probability of breaking for a fraction of defects 
P < Pc,oo and a breaking probability equal to 1 for p > Pc,oo^^- However, for networks of finite 
size there is a non vanishing probabihty of formation of the percolating cluster of defects 
(thus, of breaking) also when the defect fraction is p < Pc,oo and, for contrast, a probability 
less than 1 for p > Pc,oo^^- For this reason, when networks of finite size are considered, Pc is 
defined as the average over the statistical ensemble of the minimum values of p corresponding 
to the formation of at least one percolating cluster^^, as anticipated in Sec. II. In the case 
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oi N X N networks with N finite, it has been found^^ that (pc — Pc,oo) — cN"^/" , where v 
is the correlation length exponent (with universal value i/ = 4/3 in two-dimensions), while 
the proportionality constant, c, depends on the lattice structure. In particular, for networks 
with square- lattice structure, it has been found^^ that c ~ and thus Pc ~ Pc,oo = 0.5, 
independently of A^. In the case of rectangular Nw x A^^ networks with square-lattice 
structure, we have found that for a fixed value of Nw, Pc decreases with iV^, by reaching its 
minimum value, pw,oo-i when — > oo. More precisely: (pc — Pw,oo) ~ ^l^^'^'^-, where pw,oo 
decreases with N^. The reduction of Pc when increasing the length, by keeping constant 
the width of the network, can be explained as a consequence of the existence of a higher 
number of possible paths of defects connecting the top and the bottom of such a network 
compared to a square network with the same width. These feature is associated with the 
greatest instability of networks with this geometry. 

We note that Eqs. (9)-(ll) and the subsequent discussion apply to the case of random 
percolation, therefore the average time to failure given by Eq. (10) is independent of the 
current. On the other hand, electromigration is a current driven phenomenon which, within 
our approach, is described by a biased percolation^^'^^. The effect of the biased percolation 
can be roughly decomposed in two components : i) a correlated growth of the defect pattern, 
which exhibits an increasing degree of filamentation at increasing currents; ii) an average 
heating of the network due to Joule effect. 

The filamentation driven by the external current implies a bias dependent percolation 
threshold Pc{I), as shown in Fig. 9. In this figure the full circles represent the values of the 
percolation threshold calculated by averaging the fraction of defects which corresponds to 
the failure of 20 networks of sizes 12 x 400 subjected to external currents / > Ib- The gray 
region in Fig. 9 corresponds to values I < Ib (stationary region), while the solid curve is 
the best-fit with a quadratic expression. We notice that for low currents the defect pattern 
is found to exhibit a relatively week degree of filamentation and the value of Pc — 0.21 is 
not very far from the random percolation threshold (which, for these values of Nw and Nl, 
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is Pc — 0.37). Then, for currents up to / ~ 40 mA the filamentantion achieves its maximum 
level and consequently Pc its minimum value. At further increasing / values, the onset of a 
multi-filament ation pattern is observed^^ with the simultaneous growth of several filaments 
of defects (voids) elongated perpendicular to the direction of the current flow^^. This effect 
manifests itself in a smooth increase of the percolation threshold at the highest currents, as 
shown in Fig. 9. 

By introducing the above dependence of Pc on I in Eqs. (10)- (11) we obtain the average 
time to failure versus current shown in Fig. 10 (up-triangles) . For the sake of comparison, 
we also show in this figure (open circles) the values of ^50 obtained by MC simulations and 
already reported in the inset of Fig. 8. We can see that at low currents the ATTFs obtained 
from Eq. (10) by accounting for the bias driven filamentation, tend to merge with the MC 
results. Then, in the moderate current region the ATTF exhibits a power-law behavior, 
with a current exponent n — 1.5 (the dashed curve in Fig. 10 is the best-fit with such a 
power- law). Finally, for currents I ^ AQ mA the ATTF nearly saturates. The discrepancy at 
intermediate and high currents between the two sets of data (ATTFs from Eq. (10) and MC 
simulations) can be explained in terms of Joule heating effects. These effects can be included 
in Eq. (10), in the spirit of a mean field theory, by replacing in the random percolation 
expressions of the breaking and recovery probabilities, the temperature To with an average 
temperature < T >— Tq + AT, where AT is given by Eq. (5). By introducing into Eq. (10) 
these new average probabilities < Wqp > and < Wr > and by simultaneously accounting 
for the dependence Pc{I), we obtain the average time to failure reported as down-triangles in 
Fig. 10. Thus, the up triangles in this figure correspond to a value of the structure thermal 
resistance ^ = and the down triangles to ^ 7^ 0. The excellent agreement found between 
MC results and those obtained through Eq. (10) supports the interpretation suggested here. 
Thus, the generalization of Eq. (10) to the case of biased percolation, made by accounting 
for both Joule heating effects and for the dependence of the percolation threshold on the 
bias, is able to describe quite well both the results of the MC simulations and the behavior 
observed in many EM experiments^^^'^. 
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The results reported in Fig. 10 shed new hghts on the role played by Joule heating 
effects in determining the value of the current exponent in the Black's law. Indeed, though 
many accelerated EM tests^^^ provides a value n = 2 other values of n have been fre- 
quently measured Values of n greater than 2 have been usually attributed to Joule 
heating^"^'^^. In this respect, Fig. 10 points out quite well the importance of this effect in 
the high current region, thus confirming the interpretation made in the literature. However, 
Fig. 10 shows that Joule heating can affect the value of n also in the moderate current 
region. In this region values 1 < n < 2 have been reported by different authors^'^'^^. Many 
EM models have been proposed to explain the value of the current exponent^'^'^'''^'^. These 
models fall into two main categories: "void growth" and "nucleation" models^^. In the void 
growth models, the failure is taken to occur after a void grows up to a critical size. It is 
generally accepted that this category of models provides a value n = 1^^. In nucleation 
models, the failure arises from the buildup of a critical vacancy concentration. It is gen- 
erally beheved that this second category of models provides a value n — 2^"^. Recently, 
Tammaro and Setlik^^ have proved that even when nucleation is the limiting process it can 
be 1 < n < 2. Our results confirm this conclusion. As a matter of fact that the present 
model, which relates the failure to the achievement of a percolation threshold for the defect 
fraction, belong to the second category of models. In fact, as shown by Fig. 10, metallic lines 
with different structure thermal resistance 9, thus exploiting different sensitivity to Joule 
heating effects can exhibit different current exponent values, even in the moderate current 
region. Concerning this point, we underline the fact that a power-law behavior of the ATTF 
versus current is predicted by the model merely as a consequence of the filamentation of 
the damage pattern, even neglecting the effect of the average heating of the metallic line. 
Moreover, the dependence of n on the length of the lines, reported in EM tests^, can be ex- 
plained in terms of dependence of both, heating effects and Pc, on the length of the system. 
Finally, from Fig. 10 we can see that Joule heating is also responsible for the shift towards 
lower current values of the crossover between the high-current and the moderate-current 
regions (in this case from 40 mA to 30 mA) . On the basis of the above considerations we can 
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now understand the dependence with current of the shape factor of the TTFs distribution 
shown in Fig. 7. At relatively low currents, just above the Ib threshold, the degree of 
filamentation of the defect pattern is week and there is a wide spread in the network micro- 
geometries and thus in the TTFs. Then, at increasing currents, the degree of filamentation 
increases and the effect of the different network microgeometries is reduced together with the 
spread in the TTFs. Finally, in the highest current region, because of strong Joule heating, 
multi-filamentation occurs^^. This implies a high degree of stochasticity and in turn a high 
variabihty of microgeometries together with an increase of the TTFs spreading. 

C. Geometrical effects 

In this section we investigate the effect of the length and width of the network on the 
failure process. First, we have analysed the dependence of ^50 on the length A^^ of the 
network. Figure 11 displays the simulated median time to failure as a function of Ni^. The 
data indicated by full circles are obtained by taking: Nw — 12 and / = 10.8 mA, while the 
data reported as open squares and shown in the inset of Fig. 11 correspond to Nw — 36 
and / = 32.4 mA. Thus, in both cases the current density (measured in current units)^^ is 
j — I/Nw — 0.9 mA. We have found that t^o sharply increases by decreasing Nl and it 
diverges for network lengths below a certain value, N^^. This length can be considered as 
a critical length of the network and it is dependent on the network width, as shown in Fig. 
11. Furthermore, for sufficiently long networks, ^50 nearly saturates to a value independent 
of the length and increasing when increasing the width. In the following, this asymptotic 
value of the median time to failure in the limit of infinitely long lines will be denoted as tinf. 
The behavior of ^50 shown in Fig. 11 is in qualitative good agreement with the behavior 
observed in the EM experiments^. 

To determine more precisely the dependence of the median time to failure on the geomet- 
rical parameters, we have analysed the simulated values of t^o as a function of the difference 
{Nl — NlJ. Figure 12 reports the log- log plot of the difference t^o — tin/ as a function of the 
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difference {Nl — N^^y Here, together with the data aheady shown in Fig. 11 we have also 
reported the data for Nw = 48 and / = 43.2 mA (open circles). In this way, all the data 
in Fig. 12 correspond to the same value of j. We have found the following values of critical 
length: Lc — 2.5, 7.5, 10.0, for Nw — 12, 36, 48, respectively. Therefore, in all the cases it is: 
(NlJNw) = 0.21 ±0.01. The values oiUnf range between 2 x 10^^3 x 10^. Figure 12 shows 
that the difference t^o — Unf follows a power-law as a function of the difference {Nl — Nl^); 
indeed, the solid, dashed and long dashed curves in this figure fit the data with a power-law 
with exponent —A, where A — 0.62 ± 0.02. In particular, as shown in the inset of Fig. 12, 
we have found that the three sets of data collapse onto the same curve once the difference 
^50 ~ Unf is considered as a function of the normalized quantity {Nl — NlJ/Nl^- Thus, the 
difference ^50 — Unf scales with the ratio {Nl — NlJ/Nl^ : 



rNL-NL^l 


-A 




-A 











(12) 



With the aim of identifying the physical parameters determining the value of A, we have 
investigated the dependence of ^50 on the length of the networks when these are stressed by 
different current densities. The results of simulations are reported in Figs. 13 (a) and (b). 
Figure 13(a) displays the difference ^50 — Unf versus {Nl — NlJ for Nw — 36 and j — 0.9 
mA (open squares) and j = 0.75 mA (up-triangles) . The first set of data is the same of Fig. 
12 (thus Nl,:, = 7.5), while for the second set it is Nl,^ = 8.9 and therefore NlJNw = 0.25. 
First we note that it is: 

{jNlJi ^ {jNLj2 (13) 

in agreement with the Blech's law^'^'^^ (in this case we have obtained for the threshold 
product a value 6.8 mA ). Second, we have found that the value of A is different for the two 
sets of data. Precisely, for the data obtained by taking j — 0.75 mA it is A = 1.04 ± 0.03. 

Figure 13(b) displays the difference t^o — Unf versus {Nl — Nlc) for N!^ — 48. The 
current densities j = 0.9 mA (open circles) and j = 0.75 mA (down-triangles) are the same 
of those in Fig. 13(a), while the values of the critical lengths are respectively: A^^^^ = 10.0 

24 



and A^l^^ = 12.0. The value of the exponent found for j = 0.75 mA is now A = 1.20 ± 0.04. 
This result suggests a dependence of the exponent A not only on j. We note that also for 
this set of data the Blech's condition, Eq. (13), is satisfied with {jNi)c = 9.0 mA. Thus, 
by considering the threshold product as a function of the different parameters, {jNi)c — 
F{ro, NiY: Er,Tq...), F is found to be an increasing function of Nw 

According to the model proposed by Blech^'', the threshold product is determined by 
the ratio {Q,a^cr) I {pZ*e) (Eq. (2) in Sec. I). The quantity {pZ*e) is related to the intrinsic 
properties of the material and, within our model, can be associated with the parameter ro, 
defined in Sec. II. On the other hand, the product (riaAcr), being a function of the geometry 
of the line, of the properties of the electrical contacts, of the presence of passivation layers 
and of the temperature^'^'^'^^, within our model is controlled by the geometry of the network, 
by the efficiency of the recovery process (i.e. the energy Er) and by the temperature Tq. In 
particular, according to Eq. (2), the value of the threshold product (j-L)c is an increasing 
function of Au which is the maximum value of the mechanical stress difference between the 
line terminals that the line can bear^'^'^''. On the other hand, the value of this maximum 
stress difference increases with the fine width. Therefore, the dependence of the threshold 
product on the width predicted by our model is in qualitative agreement with the behavior 
described by Eq. (2). In any case, further studies are necessary to determine the dependence 
of tinf, A and of the threshold product on tq, Nw, Er, Tq, and the other parameters of the 
model. 

Below we analyse in more detail the dependence of the simulated median time to failure 
on the width of the network. To this purpose Fig. 14 reports t^o versus for networks of 
different lengths, stressed by a current density j = 0.9 mA. Going from the top to the bottom 
of the figure the different sets of data correspond to: A''^ = 15, 18, 36, 48, 100, respectively. 
The different curves in Fig. 14 fit the corresponding data with the expression: 



where F{x) = (1 - e"^), A = 0.62, Nl^ = 0.21 A^^;,,, Wq = 7 ^ 12, r is a constant related to 
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(14) 



the value of Unf, and X is a fitting constant. As a general trend, t^o increases systematically 
at increasing the width of the network. However, in the case of long lines the dependence of 
^50 on the width shows a saturation at the largest width. By contrast, in the case of short 
lines this saturation disappears and ^50 exhibits a final sharp increase with the width and a 
tendency to diverge when approaching the Blech condition. This behavior of ^50 is in overall 
agreement with the results of EM tests performed on lines larger than about 2 /xm^'^. On 
the other hand, it has been found that narrow lines, with width smaller than 2 /im, exhibit 
a sharp lengthening of the median time to failure This phenomenon occurs because for 
such narrow lines (known as bamboo structures) the grain size becomes comparable with 
the line width^'^. As a consequence, the lack in these lines of grain boundaries along the 
direction of current implies that EM can only occur within the bulk of the grains or along 
interfaces. These processes usually require an activation energy significantly higher than 
that associated with EM along grain boundaries''^. Therefore, bamboo structures exhibit 
median time to failure longer than that displayed by other kinds of fines'"^. Figure 14 shows 
that our model, in the present formulation, is unable to describe the lengthening of ^50 for 
very narrow lines. This is not surprising considering the fact that here we have taken a single 
value for the activation energy of the breaking process, Eqp, that is a value common to all 
the resistors in the network. A further implementation of the model which introduces two 
different activation energies for the breaking of bulk resistors and surface resistors, Eqp^b 
and Eqp^s respectively, with Eqp^b > Eqp^s, would account also for this lengthening of t^o 
for narrow lines typical of bamboo structures. 



D. Saturation and resistance fluctuations 

In this section we consider the situation occurring at low current densities or for short 
line lengths, when the product of the current density and of the line length is lower than the 
threshold value. In this case electromigration stops and the resistance of the fine achieves 
a steady state value (saturation value) dependent on the external conditions and on the 
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properties of the line^^'^. The steady state is characterized by fluctuations of the resistance 
around the average value < R > (this value is calculated by averaging over all the steady 
state values of the resistance, i.e. the values taken after the transient time associated with 
the termination of the EM process)^^. The study of resistance saturation effects and of 
their dependence on the stress current, temperature, geometry and other properties of the 
metallic line, can provide an important tool to investigate EM phenomena^, alternative to 
the study of the median time to failure. Actually, studies of EM based on saturation effects 
exhibit the advantage to be non-destructive and in particular, they reveal their effectiveness 
in the analysis of short lines^. Furthermore, from a fundamental point of view, they provide 
the possibility to investigate fluctuation phenomena under far from equilibrium conditions. 

Below we consider the steady state of networks biased by currents below the electrical 
breakdown threshold I < Ib- Figure 15 shows the resistance evolution of a 12 x 400 network 
with the same parameters considered in Sec. III.B. More precisely, we report the difference 
R{t) — Rq versus time for a network stressed at Tq — 492 K by the current 7 = 4.0 mA which 
corresponds to the threshold for electrical failure. Saturation to < i? > and large resistance 
fluctuations are well evident in Fig. 15 for this current value. The inset in this figure reports 
on an enlarged time scale the initial transient values oi R — Rq, evidencing also the initial 
decrease of the resistance due to alloying effects, discussed in Sec. II. Details concerning the 
dependence of < > on the value of the stress current, on the bias conditions (constant 
current or constant voltage) and on the TCR can be found in Refs. 49. 

Here we discuss two important features of the resistance fluctuations occurring in a 
metallic line in a non-equilibrium steady state, stressed by a current near the EM threshold. 
First, we consider the non-Gaussianity property^^~^° of the distribution oiSR = R — < R>. 
To this purpose, we have calculated the probability density function, $, of the distribution of 
5R for the steady state signal in Fig. 15. Figure 16 reports the product $(7 as a function of 
{< R > —R)/a in a lin-log plot, where a is the root mean square deviation from the average 
resistance. For comparison, in the same flgure we also report the Gaussian distribution 
(dashed curve), which in this normalized representation has zero mean and unit variance. 
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This representation has been adopted because, by making the distribution independent of its 
first and second moments, it is particularly convenient for exploring the functional form of a 
distribution®^. The results in Fig. 16 show a considerable non-Gaussianity of the distribution 
of 5R, which is associated with the fact that the system is close to breakdown conditions. 
We have found that this non-Gaussianity is sufficiently well fitted by a generahzed form®'^ 
of the Gumbel distribution (solid curve), where this last is commonly used for the analysis 
of extreme events^. The role of the stress current and of the size of the system on the 
distribution of the resistance fluctuations has been studied in Refs.®^, where the conditions 
under which the distribution achieves the universal behavior described by the Bramwell, 
Holdsworth and Pinton distribution®'^ have been identifled. 

As a second feature, we consider the power spectral density of resistance fluctuations. 
To this purpose. Fig. 17 displays the power spectrum associated with the steady state 
signal in Fig. 15. Both the frequency and the spectral density are expressed in arbitrary 
units. Two regions can be identifled in the spectrum: a 1/f like branch at low frequencies 
and a Lorentzian cut-off at high frequencies. This feature is here interpreted as due to 
the presence of three relaxation mechanisms: i) the slowest one related to the breaking and 
recovery of the regular resistors and describing the generation and recovery of microvoids; ii) 
the fastest one associated with the transformation of regular resistors into impurity resistors 
(and vice versa) and describing alloying effects; iii) an intermediate one corresponding to 
generation and breaking of impurity resistors. Accordingly, the Lorentzian branch describes 
the evolution occurring on a fast scale when only the second transformation is present. 
On the slow scale, the other two transformations become of importance and thus all the 
relaxation mechanisms coexist. Therefore, a 1/f-like spectrum appears in the low frequency 
region, in agreement with resistance noise measurements performed on Al-Cu alloys^'®'^^. We 
notice that by taking different activation energies for bulk and surface resistors, as described 
at the end of Sect. III.C, a further relaxation time would be added. 
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CONCLUSIONS 



We have applied the biased percolation model to the study of degradation and failure 
phenomena induced by EM in metallic line. Our "coarse grain" approach focuses on the 
correlations established by the electronic current among the different elemental resistors of 
the network which mimic the structural components of the system (grains, clusters of grains, 
interfaces, etc. i.e. atomic transport channels). This approach provides a unified theoret- 
ical framework able to account successfully for many relevant features of the experiments, 
including the damage pattern, the resistance evolution, alloying effects and the statistical 
properties of TTFs. In particular, the model correctly predicts a lognormal distribution 
for TTFs, perfectly superimposing with the experimental one and it is able to estimate the 
dependence of the shape factor of the distribution on current and temperature. In what 
concerns the dependence of the median time to failure on the stress conditions, the model 
predictions agree with the experiments over the full range of current and temperature values 
considered. Simulations performed on rectangular networks of different length and width 
have allowed us to investigate the dependence of TTFs on these parameters. The results of 
the model agree with the existence of a Blech's length, moreover they predict the existence 
of a scahng relation between the MTF and the hne length. Finally, we have considered 
resistance saturation effects. In this case we have studied the properties of the resistance 
fluctuations, by focusing on the non-Gaussianity of the distribution and on their power 
spectrum. The flexibility of the theoretical approach offers further possibility to describe 
and interpret phenomena which at present have not been considered. We finally emphasize 
that EM, which occurs in granular materials, in presence of a significant disorder, driven 
by an external bias and contrasted by growth of mechanical stress gradient, represents a 
paradigmatic example of failure process in a disordered system. Therefore, the ability of 
our approach to account for a wide scenario of the EM related phenomenology should be 
of interest in the more extensive perspective of understanding non-equilibrium and failure 
phenomena in disordered systems. 
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FIGURE CAPTIONS 



Fig. 1 (a). Typical resistance evolution of a 12 x 48 network. The values of all the other 
parameters are specified at the end of Sec. II. The stress conditions are: / = 10.8 mA and 
To — 492 K. The resistance is expressed in Ohm and the time in arbitrary units corresponding 
to the number of iteration steps, (b) Damage pattern at the iteration step t — 2150 of the 
evolution shown in (a). The different gray levels, from black to white, are associated with 
different T„ values, ranging from 492 to 700 K with a step of 10 K. 

Fig. 2. Experimental resistance evolutions of seven Al-0.5%0.5Cu lines stressed at T = 467 
K by a current density j = 3 MA/cm^ (which corresponds to / = 10.8 mA within the 
model). The resistance is expressed in Ohm and the time in seconds. 

Fig. 3. Calculated resistance evolutions of the Al-0.5%0.5Cu lines in Fig. 2. The simulations 
have been performed by taking Tq = 467 K and / = 10.8 mA. The value of the remaining 
parameters are specified in the text. The resistance is expressed in Ohm and the time in 
seconds by using the value At = 185 s for the time interval associated with each iterative 
step (see text) . 

Fig. 4. Lognormal plot of the cumulative distribution function of the failure probability 
(expressed in percentage) as a function of the time to failure: (full circles) TTFs experi- 
mentally measured and (open circles) calculated by the model. The data correspond to the 
same statistical samples considered in Fig. 2 and 3, respectively. The stress conditions are 
To = 467 K and / = 10.8 mA. The dashed line fits the CDFs with a lognormal distribution. 

Fig. 5. Lognormal plot of the cumulative distribution functions of the failure probability 
(expressed in percentage) as a function of the time to failure. The different functions are 
calculated at different substrate temperatures: Tq = 800 K (triangles left). To = 650 K 
(open squares), Tq = 467 K (open circles), Tq = 400 K (triangles down). The stress current 
is / = 10.8 mA. The solid lines fit the CDFs with lognormal distributions. 
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Fig. 6. Median time to failure, ^50, as a function of the inverse substrate temperature. 
The median times to failure are expressed in arbitrary units and the temperature in K. 
The stress current is / = 10.8 mA. The dashed line is the fit with the exponential function 
Z exp[4700/ro]- 

Fig. 7. Lognormal plot of the cumulative distribution functions of the failure probability 
(expressed in percentage) as a function of the time to failure. The different functions are 
calculated at different stress currents: / — 7.5 mA (triangles down), / — 30.0 mA (open 
squares), / = 60.0 mA (open circles). The substrate temperature is Tq = 492 K. The solid 
lines fit the CDFs with lognormal distributions. 

Fig. 8. Log-log plot of tso versus I — lb-, where Ii, is the breakdown current defined in the text. 
The two lines of slope —1.5 and —5.2 represent the fits with a power- law in the moderate 
current and in the high current regions, respectively. The inset shows the log-log plot of t^Q 
versus /. In both the main figure and the inset, the median times to failure are expressed 
in arbitrary units and the current in mA. The gray region in the inset corresponds to the 
stationary region attainable for currents lower than the Ib value. 

Fig. 9. Percolation threshold for broken resistors, Pc, versus current (this last is expressed 
in mA). The curve is a quadratic fit (see text). The gray region evidences the stationary 
region. 

Fig. 10. Log-log plot of versus /. The median times to failure are expressed in arbitrary 
units and the current in mA. The open circles (from MC simulations) represent the same 
data reported in the inset of Fig. 8. The up-triangles are obtained by Eq. (10) but taking 
into account the dependence Pc{I) shown in Fig. 9. The down-triangles are obtained with the 
same procedure but replacing the probabilities Wdo and Wrq with < Wqp > and < Wr >. 
The solid, long-dashed and dotted curves represent the best-fit with a power-law with slopes 
—2.1, —1.5, —5.7, respectively. 
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Fig. 11. Median time to failure, tso, as a function of the network length Nl for Nw — 12 
and Nw — 36 (inset). The median time to failures are expressed in arbitrary units. The 
dashed lines are a guide to the eyes. 

Fig. 12. Log- log plot of the difference t^o — tin/ versus the difference {N^ — Nl^^, where 
tinf is the median time to failure in the limit of infinitely long lines and Nl^ is the critical 
length. Both quantities are expressed in arbitrary units. The full circles are obtained by 
taking — 12, I — 10.8 mA. In this case is Nl^ — 2.5. The open squares correspond to 
Nw = 36, 7 = 32.4 mA and Nl, = 7.5 while the open circles to Nw = 48, 7 = 43.2 mA 
and Nl, = 10.0. The sohd, dashed and long dashed curves fit the data with a power-law of 
exponent —0.62 ± 0.02. In the inset the same data of ^50 — Unf are reported as a function of 
{Nl - Nl,)/Nl,. 

Fig. 13 (a). Log- log plot of the difference ^50 — Unf versus the difference {Nl — Nl,), where 
tinf is the median time to failure in the limit of infinitely long lines and A^^^ is the critical 
length. Both quantities are expressed in arbitrary units. The network width is A''^^ = 36. 
The open squares are obtained by taking I — 32.4 mA, in this case Nl, — 7.5. The up- 
triangles correspond to J = 27.0 mA and Nl, = 8.9. (b) Nw = 48; open circles: / = 43.2 
mA and Nl, — 10.0; down-triangles: / — 36.0 mA and Nl, — 12.0. 

Fig. 14. Plot of as a function of network width Nw- The different sets of data correspond 
to networks of different length: A^^ — 15 (down triangles), Nl — 18 (open squares), Nl — 36 
(full squares), Nl = 48 (up triangles), Nl = 100 (stars). The curves fit the corresponding 
data with the expression specified in the text. 

Fig. 15. R{t) — Rq versus time for a network 12 x 400 stressed at Tq = 300 K by a current 
I = Ib = 4: mA. Here Rq is the perfect network resistance. The inset shows on an enlarged 
time scale the evolution of the resistance in the initial stage. 
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Fig. 16. Normalized probability density function, F = ^a, of the resistance fluctuations 
reported in Fig. 16. The sohd curve fits the data with a generalized Gumbel distribution 
(see the text), while the dashed curve is the Gaussian distribution. 

Fig. 17. Power spectral density of the resistance fluctuations in Fig. 16. The frequency and 
the spectral density are expressed in arbitrary units. 
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